On the existence and scaling of structure functions in turbulence according to the data 

Alexandre Arenas 

Departament d'Enginyeria Informatica i Matematiques 
Universitat Rovira i Virgili, 43007 Tarragona, Spain 
and 

Alexandre J. Chorin 

Department of Mathematics 
University of California and Lawrence Berkeley National Laboratory 
Berkeley, CA 94720 

Abstract 

We sample a velocity field that has an inertial spectrum and a skewness that matches experimental 
data. In particular, we compute a self-consistent correction to the Kolmogorov exponent and find 
that for our model it is zero. We find that the higher order structure functions diverge for orders 
larger than a certain threshold, as theorized in some recent work. The significance of the results for 
the statistical theory of homogeneous turbulence is reviewed. 

1 Introduction 

In 1941 Kolmogorov 1 formulated his famous scaling theory of the inertial range in turbulence, according 
to which the second order structure function, i.e., the function ^(r) =< ((u(x + r) — ?i(x)) 2 >, scales as 
(er) 2 / 3 , where x, x+r are points in a turbulent flow field, u is the component of the velocity in the direction 
of r, e is the mean rate of energy dissipation, r is the length |r| of r, and the brackets denote an average. 
A Fourier transform yields the Kolmogorov-Obukhov inertial range spectrum E(k) = Ce 2 / 3 fc~ 5 / 3 , where 
C is a constant and k is the wave number |2j. These results are the lynchpins of turbulence theory yet 
uncertainty lingers as to their general validity and as to the details of the derivation. Subsequently ( 
see e.g. [3]) it was claimed that this scaling result generalizes to structure functions of any order, i.e., 
to S p (r) =< ((it(x + r) — u(x)) p >=(e7-) p / 3 for any p > 0, but note that Kolmogorov and Obukhov 
themselves never went that far; Kolmogorov in 0] gave only the further result for p = 3. 

Shortly after the publication of this theory, Landau [5] challenged its derivation on the ground that 
the rate of energy dissipation is intermittent, i.e., spatially inhomogenous, and cannot be treated as a 
constant. This observation fits in with the experimental observation that for p > 3 the exponents are 



smaller than the values given by the extended Kolmogorov theory and, and it is also often claimed for 
the p = 2 the experimental value of the Kolmogorov-Obukhov exponent is larger than the prediction of 
this scaling theory. Kolmogorov and Obukhov themselves produced a "corrected" theory [HI, Eli an( i the 
current belief is that the theory has to be supplemented by "intermittency corrections" ; much effort has 
been expended on the calculation of these corrections (see e.g. 

A different analysis of possible corrections to the Kolmogorov exponents has been offered by Barenblatt 
and Chorin in |SJ . The structure functions depend on the Reynolds number R (for a definition of 
R suitable for the present case see a bulk length scale L, and the mean rate of energy dissipation 

e. Dimensional analysis yields S p = (er) p / 3 $ p (r/L, R), where <& p is an unknown dimensionless function 
of two large arguments. If one makes a complete similarity assumption (see )12|). one finds & p ~ C for 
large arguments, and one recovers the Kolmogorov exponents above, but other assumptions are possible 
and indeed natural. In particular, an analogy with a successful scaling theory for turbulent boundary 
layers ^Sjj |14| . leads to the assumption <E> P = C p (R){r / L)( ap / lnR \ where a p is a constant and C P (R) is a 
function of the Reynolds number R. This yields structure functions S p proportional to r p/ 3 + a p/ tnR : j. e ., 
exponents that depend on R, and maybe even more important, a proportionality constant that depends 
on R and may well diverge as R — -> oo for p large enough. As pointed out in the dependence on 
1/lnR is an assumption, and in reality may be weaker yet; such dependence could be hard to detect in 
experimental data. 

This analysis is quite compatible with the Landau observation regarding intermittency, as was al- 
ready pointed out in |14j : Turbulent flow is intrinsically intermittent, and indeed in the limit of infinite 
R, extremely intermittent, in the sense that the bulk of the vorticity occupies a negligible fraction of 
the available volume |15|: this is a dominant fact about turbulence, and the notion of "intermittency 
corrections" is no more meaningful than a notion of "turbulence corrections" to turbulent flow. However, 
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at finite R viscosity reduces the intermittency and the scaling has to be corrected for "intermittency 
reductions" that depend on R. In related work |16| . Barenblatt and Chorin made numerical estimates 
on the basis of an approximate theory of turbulence and came to the conclusion that, as R increases, the 
structure functions for p > 4 diverge; similar conclusions had been reached by P. Constantin (personal 
communication) through mathematical considerations and by Mandelbrot |17| by a geometrical similarity 
argument. Thus indeed C P (R) would diverge as R increases for p > 4. 

In the absence of solutions of the Navier-Stokes equations these various theories can only be checked 
by data from experiments, and this seems to be difficult at present. It occurred to us to interrogate 
further the available data about the inertial range by building a stochastic computer model of the inertial 
range, compatible with well-accepted and reproducible data about their skewness (indicating the level of 
intermittency), and then examine the resulting structure functions. The field we construct is obviously 
not unique; however, any model that satisfies constraints imposed by the data is highly instructive, and 
indicates what possibilities are worthy of further study. The results are striking enough to open new 
vistas for theory, as we discuss in the concluding section. 

We thus construct on the computer a one-dimensional Gaussian velocity field with a power-law, then 
modify it so that the velocity differences assume the skewnesses presented in Batchelor's book ^S] on the 
basis of the data in [T3\ ; the power law will be either the Kolmogorov-Obukhov spectrum or a modification 
that satisfies other constraints (see below). We use a one-dimensional model so that we can perform the 
calculations with sufficent accuracy, and think that this model is sufficient to carry our conclusions. We 
then calculate, when possible, the higher-order structure functions. The data in |18l IH)] are tabulated 
there for various Reynolds numbers; we take the data that correspond to the largest R = 42200; our model 
has a power law for all wave numbers and is thus inviscid. Note that by modern standards this value of 
R is not very large; we thus assume that the skewness does not vary drastically as R increases further. 
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We first present the technical aspects of the construction, then the results, then provide a discussion. 

2 Building a Model 

We first explain how to sample effectively a homogenous Gaussian velocity field with a given spectrum. 
Our basic tool is a construction due to Elliott et al. |2U1[^. A homogeneous Gaussian random field is 
determined by its second order means: 



< u(x) > = 0, 

< u(x + Ax)u(x) >= I e- 2mkx E{k)dk. (1) 



The spectral representation of the Gaussian random field u is 



u(a;) = / e 2mkx E 1 / 2 (k)dw(k), (2) 

J — OO 

where it; is a a Wiener process. Expand dw in a complete orthonormal series </> m : 

OO 

dw(k) = 

(k)dk, (3) 

m— 1 

where the 7 m form a sequence of independent random Gaussian variables. By substituting © into Eq. 
© we find: 

u{x) = 7mC m (l), (4) 

m 

where the coefficients c m (x) are: 

Cm(x)= e 2mkx E 1 ' 2 {k)4> m {k)dk = T- 1 [E 1 ' 2 <t> m }(x), (5) 
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and T is the Fourier transform. Elliott et al. proposed to use as basis of the decomposition, Fourier 
transforms of wavelets based on the Meyer wavelet; these wavelets are generated from the Meyer mother 
wavelet ip(k) by the wavelet relation: 

4> mn (x) = 2- m ^(2- m x-n), (6) 

where the index m refers to different scales (octaves) and n to different dilations. Using this prescription, 
Eq. (Jijl can be represented as, 



m n 

c mn (x) -^- 1 [</ 2 0](2-™x-n), (7) 

where ££/ 2 (fc) = 2- m / 2 E 1 / 2 (2-' m k). In the particular case of the Kolmogorov-Obukhov energy spectrum 
E{k) ~ fc~ 5 / 3 , the Fourier transform of the c mn coefficients represented in Eq. Q becomes: 

(Fc mn )(k) = 2-" l / 2 | 2-™/ 2 fc f 5/ V mn (fc). (8) 
Similar formulas are obtained for every power-law spectrum. 

Defining f m {2' m x - n) = F- 1 (2- m / 2 \ 2- m / 2 k f 5/6 ?/wO)), the representation of the field given in 
Eq.© is: 

< x ) = EE 7™.J m (2" m ^ - n). (9) 
The use of the summation (0 requires a truncation in m and n that preserves accuracy. The truncation 
uses the spatial decay of the functions f m (see \2Q\). It is convenient to center the summation around the 
term with the smallest (in magnitude) value of the argument (2~ m x — n). This can be done by defining 



for each term m of the outer sum the index n m = \2~ m x\ and shifting indexes in the inner sum to 
n' = n — h m . The parameter m governs the scale, the truncation in m defines a frequency range, which 
can be modified by rescaling the distance x if necessary, so as to obtain a numerically convenient range 
— M < m < 0. The truncation in n (the number of dilations) defines the region where the support of f m 
is concentrated; we suppose that this range lies between a pair of integers -N and N . We can then write: 

N 

U ( X ) = E E lrnA m+ nfm{2~ m X-[2- m x\-n). (10) 

m=— M n=—N+l 

More information can be found in (201E21 A technical issue concerning the generation of Gaussian random 
variables should be pointed out. The number of random variables needed to sample the field scales as 
2 M+1 x, and this imposes severe demands on the storage in the computer. To overcome this problem, 
a reversible pseudorandom number generation is necessary to sample "f mn . Simple linear congruential 
generators with moduli around 2 31 can exhaust their period in few minutes in a conventional PC, and the 
resulting poor distribution of the samples can dramatically bias simulation results for sample sizes much 
smaller than the period length. To overcome this problem, we used reversible multiple linear congruential 
generators with many long streams and substreams \2'6\ that provide periods of approximately 2 191 . 

We now explain how to modify the construction of the preceding section so as to impose on the 
sampled velocity field the non-Gaussian characteristics observed in the experimental data. 

We introduce non-Gaussianity into the numerical experiment via the coefficients 7 m . In the previous 
paragraphs the 7^s were Gaussian variables; now we construct a new distribution for the "f' m s that 
preserves the mean and variance of the original Gaussian distribution but with a skewness different from 
0. The skewness of a randon variable 77 is < if > / < rj 2 > 3 / 2 ; it is zero for a Gaussian variable. We 
describe a transformation that maps a Gaussian variable on a non-Gaussian variable with the same mean 
and variance but with a prescribed, non-zero skewness controlled by a single parameter. A Gaussian 
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Figure 1: Non-Gaussian probability density functions, obtained by the algorithm in the text 



variable with mean zero and unit variance has the probability density 



1 



p{x) 



-x 2 /2 



'2tt 



(11) 



The following change of variables yields a new skewed probability density function, with a negative 
skewness as in the cases we consider: 



2/0) 



-(e 



-ax g,a /2 



(12) 



The new probability density function g(y) is obtained by calculating g{y) = p(x)/\ provided y{x) 
is a monotonic function. This formula provides a new distribution with zero mean, unit variance and a 
skewness controlled by the parameter a. For the Ea. l|12Jl the skewness coefficient C3 reads: 



C3 



A plot of these transformed distributions g(y) is presented in Figure 1. 



(13) 
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Figure 2: Skewness as a function of separate r/M, from Stewart (^13); M is a reference length in that 
paper. Reproduced with permission. 
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Figure 3: Exponent £3 of the structure function of order 3 versus a in the relation E(k) <~ k^ +a . The 
intermittency correction should be where £3=1. 

Experimental results for the velocity field were obtained by Stewart in 1951 |19j . and are also presented 
in ^8], see Figure 2. The skewness varies with scale, and this is easy to incorporate into our wavelet 
representation as the various wavelets describe motion on different scales. The parameter a(r) has been 
obtained by inversion of the formula l|13|) for values of Ca(r) corresponding to the experiment for the 
highest value of R reported in |T5) . 

3 Results 

Having constructed these fields, i.e., written a computer program that samples them, we proceeded to 
calculate their properties by Monte Carlo; we now list some of the results. 

1. First, we noticed that if the second order structure has the Kolmogorov form, and if the field 
is Gaussian, then all the structure functions obey the Kolmogorov scaling. This is quite natural and 
obvious, and was well known to Kolmogorov (personal communication by G.I. Barcnblatt), but we had 
never seen it stated in print. The obvious converses are not necessarily true. 



2. It is "almost" a theorem that the third order structure function scales like r (see |I] as well as|15l!5]: 
all one has to assume is an extremely likely bound on the rate of blow-up of the dissipation as the Reynolds 
number grows. It is natural to ask then what form must the second order structure function take if the 
third order scales as r and the skewness is as observed. Let the second order structure function scale 
as r 2 / 3+Q and the third order structure scale as r^ 3 ; in Figure 3 we plot the computed values of £3 as 
a function of a. The relationship between them can be approximated by £3 = 1 — |a. The value of a 
should be the one that yields £3 = 1; thus, within the errors in our sampling scheme, one can construct a 
velocity field with the second order structure function with the Kolmogorov exponent, the correct third 
order exponent, and a realistic intermittency; This shows that a correction to the Kolmogorov-Obukhov 
exponent is not imposed simply by the presence of intermittency (although we cannot, of course, exclude 
that it is imposed by other constraints). Note that this analysis simply confirms that the analysis in 
Kolmogorov's paper 0] remains valid in the presence of a realistic skewness. 

3. The next question is maybe the most significant: Suppose the second order structure function has 
the Kolmogorov-Obukhov form and the skewness is as observed, what can one say about the higher order 
structure functions? When we tried to sample these higher-order functions as above, we found that the 
variance of the Monte-Carlo results was very large and did not go down as rapidly as one may expect as 
the number of trials went up, leading us to suspect that some moments of the velocity field diverged. 

We therefore plotted the data for the variable x = S\(r) = \ Au\ at fixed values of r in log-log scale (see 
Figure 4). The data for the non-Gaussian field present a power-law tail that clearly contrasts with the 
exponential decay of the Gaussian field. To estimate the form of the decay we construct the cumulative 
density function (CDF) P(x > X) of the the variable x = S±(r). We use CDF's rather than the density 
functions themselves because they fluctuate less. The behaviour of the CDFs reveals the convergence, 
or lack thereof, of the structure function moments. In the limit as the number of samples n — ► 00, 
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Figure 4: Cumulative distribution function of Si = |Au| for a Gaussian and our non-Gaussian fields at 
r /M = 1. M is again a reference length in the tabulated data. 

< S p > r — J{\ Au|) p /(| Au|)d(| Au\), where / is the probability density function, the derivative of the 
CDF. The CDF for the Gaussian field presents an exponential decay compatible with the finiteness of 
moments of any order. However the CDF for the velocity differences in the non-Gaussian case presents a 
power-law tail whose exponent determines the order of the moments that are finite. The tail of the CDF 
of Si scales for small values as a power law with exponent — 2.2±0.01, however, the final part of the CDF, 
which contain the relevant information about the decay, is well described by an exponent —3.3 ± 0.05 
within our experimental error. This means that the PDF will decay approximately as a power law with 
exponent —4.3 ± 0.05 and then produces a threshold for the higher order structure functions to converge 
at order p 3.3. We therefore conclude that, within our model, the structure functions of order > 3.3 do 
not exist. In Figure 5 we plotted the CDF's for several values of r, to show that the behavior we describe 
is at least approximately independent of r when r is small. 

This does not of course show that the moments of the true velocity field in turbulence fail to exist, 
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Figure 5: Cumulative distribution functions of Si obtained from a Gaussian field (Top), and a non- 
Gaussian field, at several separations. 

but only that the data contain a strong enough deviation from Gaussianity for these moments not to 
exist. 

Finally, we would like to provide a short discussion of sampling errors. All the data presented have 
been obtained by the least squares method, except the data for the final part of the tail in Figures 4 
and 5, where we have also applied the maximum likelihood method to reduce the estimation error of the 
exponent j23, within a framework suggested in [23|. With 95% confidence the errors are comparable 
with the thickness of the lines drawn in all cases. The regions of the sampling where aliasing problems 
appear in the FFT have been discarded in the analysis. 

4 Conclusions 

Despite the non-uniqueness of the fields we have constructed, some conclusions can be drawn from the 
computations above; in particular: 

(1) The original Kolmogorov-Obukhov —5/3 spectrum is consistent with the "exact" scaling of the 
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third order structure function even in the presence of intermittency; intermittency does not necessarily 
require that the spectrum be modified. As long as one considers only the lower order structure functions 
originally considered by Kolmogorov and Obukhov, our model provides no reason to modify their original 
scaling. 

(2) The Kolmogorov scaling of the structure functions and its extensions is exact at all orders for 
Gaussian fields; however, the data for real flows reveal enough departure from Gaussianity for structure 
functions of order higher than a threshold larger than three to diverge as the Reynolds number grows; 
if this happens, the higher-order exponents become strongly dependent on the Reynolds number, their 
scaling laws depend on the Reynolds number, and the proposals of Barenblatt, Chorin, and Prostokishin 
become eminently reasonable, as well as those of Mandelbrot and Constantin. Thus, within the limitations 
of our model, the complete similarity assumption, on which the extension of the Kolmogorov scaling rests, 
fails above the threshold but not below it; when it fails it should be replaced by an incomplete similarity 
assumption as in the papers quoted above. 
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